home *** CD-ROM | disk | FTP | other *** search
Wrap
# Source Generated with Decompyle++ # File: in.pyc (Python 2.4) import ranlib import Numeric import LinearAlgebra import sys import math from types import * class ArgumentError(Exception): pass def seed(x = 0, y = 0): '''seed(x, y), set the seed using the integers x, y; Set a random one from clock if y == 0 ''' if type(x) != IntType or type(y) != IntType: raise ArgumentError, 'seed requires integer arguments.' if y == 0: import time t = time.time() ndigits = int(math.log10(t)) base = 10 ** (ndigits / 2) x = int(t / base) y = 1 + int(t % base) ranlib.set_seeds(x, y) seed() def get_seed(): '''Return the current seed pair''' return ranlib.get_seeds() def _build_random_array(fun, args, shape = []): if isinstance(shape, IntType): shape = [ shape] if len(shape) != 0: n = Numeric.multiply.reduce(shape) s = apply(fun, args + (n,)) s.shape = shape return s else: n = 1 s = apply(fun, args + (n,)) return s[0] def random(shape = []): '''random(n) or random([n, m, ...]) returns array of random numbers''' return _build_random_array(ranlib.sample, (), shape) def uniform(minimum, maximum, shape = []): '''uniform(minimum, maximum, shape=[]) returns array of given shape of random reals in given range''' return minimum + (maximum - minimum) * random(shape) def randint(minimum, maximum = None, shape = []): '''randint(min, max, shape=[]) = random integers >=min, < max If max not given, random integers >= 0, <min''' if not isinstance(minimum, IntType): raise ArgumentError, 'randint requires first argument integer' if maximum is None: maximum = minimum minimum = 0 if not isinstance(maximum, IntType): raise ArgumentError, 'randint requires second argument integer' a = (maximum - minimum) * random(shape) if isinstance(a, Numeric.ArrayType): return minimum + a.astype(Numeric.Int) else: return minimum + int(a) def random_integers(maximum, minimum = 1, shape = []): '''random_integers(max, min=1, shape=[]) = random integers in range min-max inclusive''' return randint(minimum, maximum + 1, shape) def permutation(n): '''permutation(n) = a permutation of indices range(n)''' return Numeric.argsort(random(n)) def standard_normal(shape = []): '''standard_normal(n) or standard_normal([n, m, ...]) returns array of random numbers normally distributed with mean 0 and standard deviation 1''' return _build_random_array(ranlib.standard_normal, (), shape) def normal(mean, std, shape = []): '''normal(mean, std, n) or normal(mean, std, [n, m, ...]) returns array of random numbers randomly distributed with specified mean and standard deviation''' s = standard_normal(shape) return s * std + mean def multivariate_normal(mean, cov, shape = []): '''multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...]) returns an array containing multivariate normally distributed random numbers with specified mean and covariance. mean must be a 1 dimensional array. cov must be a square two dimensional array with the same number of rows and columns as mean has elements. The first form returns a single 1-D array containing a multivariate normal. The second form returns an array of shape (m, n, ..., cov.shape[0]). In this case, output[i,j,...,:] is a 1-D array containing a multivariate normal.''' mean = Numeric.array(mean) cov = Numeric.array(cov) if len(mean.shape) != 1: raise ArgumentError, 'mean must be 1 dimensional.' if len(cov.shape) != 2 or cov.shape[0] != cov.shape[1]: raise ArgumentError, 'cov must be 2 dimensional and square.' if mean.shape[0] != cov.shape[0]: raise ArgumentError, 'mean and cov must have same length.' if isinstance(shape, IntType): shape = [ shape] final_shape = list(shape[:]) final_shape.append(mean.shape[0]) x = ranlib.standard_normal(Numeric.multiply.reduce(final_shape)) x.shape = (Numeric.multiply.reduce(final_shape[0:len(final_shape) - 1]), mean.shape[0]) (u, s, v) = LinearAlgebra.singular_value_decomposition(cov) x = Numeric.matrixmultiply(x * Numeric.sqrt(s), v) Numeric.add(mean, x, x) x.shape = final_shape return x def exponential(mean, shape = []): '''exponential(mean, n) or exponential(mean, [n, m, ...]) returns array of random numbers exponentially distributed with specified mean''' x = random(shape) Numeric.log(x, x) Numeric.subtract(0.0, x, x) Numeric.multiply(mean, x, x) return x def beta(a, b, shape = []): '''beta(a, b) or beta(a, b, [n, m, ...]) returns array of beta distributed random numbers.''' return _build_random_array(ranlib.beta, (a, b), shape) def gamma(a, r, shape = []): '''gamma(a, r) or gamma(a, r, [n, m, ...]) returns array of gamma distributed random numbers.''' return _build_random_array(ranlib.gamma, (a, r), shape) def F(dfn, dfd, shape = []): '''F(dfn, dfd) or F(dfn, dfd, [n, m, ...]) returns array of F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator.''' return chi_square(dfn, shape) / dfn / chi_square(dfd, shape) / dfd def noncentral_F(dfn, dfd, nconc, shape = []): '''noncentral_F(dfn, dfd, nonc) or noncentral_F(dfn, dfd, nonc, [n, m, ...]) returns array of noncentral F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator, and noncentrality parameter nconc.''' return noncentral_chi_square(dfn, nconc, shape) / dfn / chi_square(dfd, shape) / dfd def chi_square(df, shape = []): '''chi_square(df) or chi_square(df, [n, m, ...]) returns array of chi squared distributed random numbers with df degrees of freedom.''' return _build_random_array(ranlib.chisquare, (df,), shape) def noncentral_chi_square(df, nconc, shape = []): '''noncentral_chi_square(df, nconc) or chi_square(df, nconc, [n, m, ...]) returns array of noncentral chi squared distributed random numbers with df degrees of freedom and noncentrality parameter.''' return _build_random_array(ranlib.noncentral_chisquare, (df, nconc), shape) def binomial(trials, p, shape = []): '''binomial(trials, p) or binomial(trials, p, [n, m, ...]) returns array of binomially distributed random integers. trials is the number of trials in the binomial distribution. p is the probability of an event in each trial of the binomial distribution.''' return _build_random_array(ranlib.binomial, (trials, p), shape) def negative_binomial(trials, p, shape = []): '''negative_binomial(trials, p) or negative_binomial(trials, p, [n, m, ...]) returns array of negative binomially distributed random integers. trials is the number of trials in the negative binomial distribution. p is the probability of an event in each trial of the negative binomial distribution.''' return _build_random_array(ranlib.negative_binomial, (trials, p), shape) def multinomial(trials, probs, shape = []): '''multinomial(trials, probs) or multinomial(trials, probs, [n, m, ...]) returns array of multinomial distributed integer vectors. trials is the number of trials in each multinomial distribution. probs is a one dimensional array. There are len(prob)+1 events. prob[i] is the probability of the i-th event, 0<=i<len(prob). The probability of event len(prob) is 1.-Numeric.sum(prob). The first form returns a single 1-D array containing one multinomially distributed vector. The second form returns an array of shape (m, n, ..., len(probs)). In this case, output[i,j,...,:] is a 1-D array containing a multinomially distributed integer 1-D array.''' probs = Numeric.array(probs) if len(probs.shape) != 1: raise ArgumentError, 'probs must be 1 dimensional.' if type(shape) == type(0): shape = [ shape] final_shape = shape[:] final_shape.append(probs.shape[0] + 1) x = ranlib.multinomial(trials, probs.astype(Numeric.Float32), Numeric.multiply.reduce(shape)) x.shape = final_shape return x def poisson(mean, shape = []): '''poisson(mean) or poisson(mean, [n, m, ...]) returns array of poisson distributed random integers with specifed mean.''' return _build_random_array(ranlib.poisson, (mean,), shape) def mean_var_test(x, type, mean, var, skew = []): n = len(x) * 1.0 x_mean = Numeric.sum(x) / n x_minus_mean = x - x_mean x_var = Numeric.sum(x_minus_mean * x_minus_mean) / (n - 1.0) print '\nAverage of ', len(x), type print '(should be about ', mean, '):', x_mean print 'Variance of those random numbers (should be about ', var, '):', x_var if skew != []: x_skew = Numeric.sum(x_minus_mean * x_minus_mean * x_minus_mean) / 9998.0 / x_var ** (3.0 / 2.0) print 'Skewness of those random numbers (should be about ', skew, '):', x_skew def test(): (x, y) = get_seed() print 'Initial seed', x, y seed(x, y) (x1, y1) = get_seed() if x1 != x or y1 != y: raise SystemExit, 'Failed seed test.' print 'First random number is', random() print 'Average of 10000 random numbers is', Numeric.sum(random(10000)) / 10000.0 x = random([ 10, 1000]) if len(x.shape) != 2 and x.shape[0] != 10 or x.shape[1] != 1000: raise SystemExit, 'random returned wrong shape' x.shape = (10000,) print 'Average of 100 by 100 random numbers is', Numeric.sum(x) / 10000.0 y = uniform(0.5, 0.59999999999999998, (1000, 10)) if len(y.shape) != 2 and y.shape[0] != 1000 or y.shape[1] != 10: raise SystemExit, 'uniform returned wrong shape' y.shape = (10000,) if Numeric.minimum.reduce(y) <= 0.5 or Numeric.maximum.reduce(y) >= 0.59999999999999998: raise SystemExit, 'uniform returned out of desired range' print 'randint(1, 10, shape=[50])' print randint(1, 10, shape = [ 50]) print 'permutation(10)', permutation(10) print 'randint(3,9)', randint(3, 9) print 'random_integers(10, shape=[20])' print random_integers(10, shape = [ 20]) s = 3.0 x = normal(2.0, s, [ 10, 1000]) if len(x.shape) != 2 and x.shape[0] != 10 or x.shape[1] != 1000: raise SystemExit, 'standard_normal returned wrong shape' x.shape = (10000,) mean_var_test(x, 'normally distributed numbers with mean 2 and variance %f' % (s ** 2,), 2, s ** 2, 0) x = exponential(3, 10000) mean_var_test(x, 'random numbers exponentially distributed with mean %f' % (s,), s, s ** 2, 2) x = multivariate_normal(Numeric.array([ 10, 20]), Numeric.array(([ 1, 2], [ 2, 4]))) print '\nA multivariate normal', x if x.shape != (2,): raise SystemExit, 'multivariate_normal returned wrong shape' x = multivariate_normal(Numeric.array([ 10, 20]), Numeric.array([ [ 1, 2], [ 2, 4]]), [ 4, 3]) print 'A 4x3x2 array containing multivariate normals' print x if x.shape != (4, 3, 2): raise SystemExit, 'multivariate_normal returned wrong shape' x = multivariate_normal(Numeric.array([ -100, 0, 100]), Numeric.array([ [ 3, 2, 1], [ 2, 2, 1], [ 1, 1, 1]]), 10000) x_mean = Numeric.sum(x) / 10000.0 print 'Average of 10000 multivariate normals with mean [-100,0,100]' print x_mean x_minus_mean = x - x_mean print 'Estimated covariance of 10000 multivariate normals with covariance [[3,2,1],[2,2,1],[1,1,1]]' print Numeric.matrixmultiply(Numeric.transpose(x_minus_mean), x_minus_mean) / 9999.0 x = beta(5.0, 10.0, 10000) mean_var_test(x, 'beta(5.,10.) random numbers', 0.33300000000000002, 0.014) x = gamma(0.01, 2.0, 10000) mean_var_test(x, 'gamma(.01,2.) random numbers', 2 * 100, 2 * 100 * 100) x = chi_square(11.0, 10000) mean_var_test(x, 'chi squared random numbers with 11 degrees of freedom', 11, 22, 2 * Numeric.sqrt(2.0 / 11.0)) x = F(5.0, 10.0, 10000) mean_var_test(x, 'F random numbers with 5 and 10 degrees of freedom', 1.25, 1.3500000000000001) x = poisson(50.0, 10000) mean_var_test(x, 'poisson random numbers with mean 50', 50, 50, 0.14000000000000001) print '\nEach element is the result of 16 binomial trials with probability 0.5:' print binomial(16, 0.5, 16) print '\nEach element is the result of 16 negative binomial trials with probability 0.5:' print negative_binomial(16, 0.5, [ 16]) print '\nEach row is the result of 16 multinomial trials with probabilities [0.1, 0.5, 0.1 0.3]:' x = multinomial(16, [ 0.10000000000000001, 0.5, 0.10000000000000001], 8) print x print 'Mean = ', Numeric.sum(x) / 8.0 if __name__ == '__main__': test()